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BACKGROUND OF THE INVENTION 

Technical Field 

The present invention relates to a system estimation 
method and program, a recording medium, and a system 
estimation device, and particularly to a system estimation 
method and program, a recording medium, and a system 
estimation device, in which the generation of robustness in 
state estimation and the optimization of a forgetting factor 
are simultaneously realized by using a fast H«> filtering 
algorithm of a hyper H w filter developed on the basis of an Hoo 
evaluation criterion . 

Background Art 
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In general, system estimation means estimating a 
parameter of a mathematical model (transfer function, impulse 
response, etc.) of an input/output relation of a system based 
on input/output data. Typical application examples include an 
echo canceller in international communication, an automatic 
equalizer in data communication, an echo canceller and sound 
field reproduction in a sound system, active noise control in 
a vehicle etc. and the like. For more information, see non- 
patent document 1: "DIGITAL SIGNAL PROCESSING HANDBOOK" 1993, 
The Institute of Electronics, Information and Communication 
Engineers, and the like. 

(Basic Principle) 

Fig. 8 shows an example of a structural view for system 
estimation (unknown system may be expressed by an IIR 
(Infinite Impulse Response) filter) . 

This system includes an unknown system 1 and an adaptive 
filter 2. The adaptive filter 2 includes an FIR digital 
filter 3 and an adaptive algorithm 4. 

Hereinafter, an example of an output error method to 
identify the unknown system 1 will be described. Here, u k 
denotes an input of the unknown system 1, d k denotes an output 
of the system, which is a desired signal, and d^ k denotes an 
output of the filter. (Incidentally, means an estimated 

value and should be placed directly above a character, 
however, it is placed at the upper right of the character for 
input convenience. The same applies hereinafter.). 

Since an impulse response is generally used as a 
parameter of an unknown system, the adaptive filter adjusts a 
coefficient of the FIR digital filter 3 by the adaptive 
algorithm so as to minimize an evaluation error e k = d k - d^ k 
of the figure. 

Besides, conventionally, a Kalman filter based on an 
update expression (Riccati equation) of an error covariance 
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matrix has been widely used for the estimation of a parameter 
(state) of a system. The details are disclosed in non-patent 
document 2: S. Haykin: Adaptive filter theory, Prentice-Hall 
(1996) and the like. 

Hereinafter, the basic principle of the Kalman filter 
will be described. 

A minimum variance estimate x A k | k of a state x k of a linear 
system expressed in a state space model as indicated by the 
following expression : 

x k+i = p~ 1/2 x k , y k = H k x k + v k (1) 

is obtained by using an error covariance matrix Z^kjk-i of the 
state as follows. 

&k\k = + K k (y k - H k x k \ k -i) 

£k+i\k = P"£fc|fc (2) 

K k = E klk -iH k T (p+ H.S^Hjy 1 (3) 

£ k \ k = E k \ k -\ — K k H k U k \ k __ Y 

£ k +\\ k = £ k \ k /p (4) 



where, 



x 0 |-i = 0, 27 0 |_i = e 0 I, £ 0 >0 (5) 



x k : State vector or simply a state; unknown and this is an 
object of estimation. 

y k : Observation signal; input of a filter and known. 
H k : Observation matrix; known. 
v k : Observation noise; unknown. 

p: Forgetting factor; generally determined by trial and error. 
K k : Filter gain; obtained from matrix Z^kjk-i- 

Z^kik: Corresponds to the covariance matrix of an error of x A k(k ; 
obtained by a Riccati equation. 

X^k+ijk: Corresponds to the covariance matrix of an error of 
x A k +i|k; obtained by the Riccati equation. 
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Z A i|o: Corresponds to the covariance matrix in an initial 
state; although originally unknown, e 0 I is used for 
convenience . 

The present inventor has already proposed a system 
identification algorithm by a fast Hoo filter (see patent 
document 1) . This is such that an Hoo evaluation criterion is 
newly determined for system identification, and a fast 
algorithm for the hyper Hoo filter based thereon is developed, 
while a fast time-varying system identification method based 
on this fast Hoo filtering algorithm is proposed. The fast Hoo 
filtering algorithm can track a time-varying system which 
changes rapidly with a computational complexity of O (N) per 
unit-time step. It matches perfectly with a fast Kalman 
filtering algorithm at the limit of the upper limit value. By 
the system identification as stated above, it is possible to 
realize the fast real-time identification and estimation of 
the time -invariant and time- varying systems . 

Incidentally, with respect to methods normally known in 
the field of the system estimation, see, for example, non- 
patent documents 2 and 3 . 

(Applied Example to Echo Canceller) 

In a long distance telephone circuit such as an 
international telephone, a four-wire circuit is used from the 
reason of signal amplification and the like. On the other 
hand, since a subscriber's circuit has a relatively short 
distance, a two-wire circuit is used. 

Fig. 9 is an explanatory view concerning a communication 
system and an echo. A hybrid transformer as shown in the 
figure is introduced at a connection part between the two-wire 
circuit and the four-wire circuit, and impedance matching is 
performed. When the impedance matching is complete, a signal 

(sound) from a speaker B reaches only a speaker A. However, 
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in general, it is difficult to realize the complete matching, 
and there occurs a phenomenon in which part of the received 
signal leaks to the four-wire circuit, and returns to the 
receiver (speaker A) after being amplified. This is an echo 
(echo) . As a transmission distance becomes long (as a delay 
time becomes long) , the influence of the echo becomes large, 
and the quality of a telephone call is remarkably deteriorated 
(in the pulse transmission, even in the case of short 
distance, the echo has a large influence on the deterioration 
of a telephone call) . 

Fig. 10 is a principle view of an echo canceller. 

Then, as shown in the figure, the echo canceller (echo 
canceller) is introduced, an impulse response of an echo path 
is successively estimated by using a received signal which can 
be directly observed and an echo, and a pseudo-echo obtained 
by using it is subtracted from the actual echo to cancel the 
echo and to remove it. 

The estimation of the impulse response of the echo path 
is performed so that the mean square error of a residual echo 
e k becomes minimum. At this time, elements to interfere with 
the estimation of the echo path are circuit noise and a signal 
(sound) from the speaker A. In general, when two speakers 
simultaneously start to speak (double talk) , the estimation of 
the impulse response is suspended. Besides, since the impulse 
response length of the hybrid transformer is about 5 0 [ms] , 
when the sampling period is made 125 [\xs] , the order of the 
impulse response of the echo path becomes actually about 400. 

Non -patent document 1 

"DIGITAL SIGNAL PROCESSING HANDBOOK" 1993 The Institute 
of Electronics, Information and Communication Engineers 
Non-patent document 2 

S. Haykin: Adaptive filter theory, Prentice-Hall (1996) 
Non-patent document 3 
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B. Hassibi, A. H. Sayed, and T. Kailath: "Indefinite- 
Quadratic Estimation and Control", SIAM (1996) 
Patent document 1 
JP-A-2002-135171 

Disclosure of tho Inv o nti on brief summary of invention 

However, in the conventional Kalman filter including the 
forgetting factor p as in the expressions (1) to (5), the 
value of the forgetting factor p must be determined by trial 
and error and a very long time has been required. Further, 
there has been no means for judging whether the determined 
value of the forgetting factor p is an optimal value. 

Besides, with respect to the error covariance matrix used 
in the Kalman filter, it is known that a quadratic form to an 
arbitrary vector, which is originally not zero, is always 
positive (hereinafter referred to as "positive definite"), 
however, in the case where calculation is performed by a 
computer at single precision, the quadratic form becomes 
negative (hereinafter referred to as "negative definite"), and 
becomes numerically unstable. Besides, since the amount of 
calculation is 0(N 2 ) (or 0(N 3 )), in the case where the 
dimension N of the state vector x k is large, the number of 
times of arithmetic operation per time step is rapidly 
increased, and it has not been suitable for a real-time 
processing . 

In view of the above, the present invention has an object 
to establish an estimation method which can theoretically 
optimally determine a forgetting factor, and to develop an 
estimation algorithm and a fast algorithm which are 
numerically stable. Besides, the invention has an object to 
provide a system estimation method which can be applied to an 
echo canceller in a communication system or a sound system, 
sound field reproduction, noise control and the like. 
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In order to solve the problem, according to the 
invention, a newly devised H*, optimization method is used to 
derive a state estimation algorithm in which a forgetting 
factor can be optimally determined. Further, instead of an 
error covariance matrix which should always have the positive 
definite, its factor matrix is updated, so that an estimation 
algorithm and a fast algorithm, which are numerically stable, 
are developed. 

According to first solving means of the invention, 

a system estimation method and program and a computer 
readable recording medium recording the program are for making 
state estimation robust and optimizing a forgetting factor p 
simultaneously in an estimation algorithm, in which 

for a state space model expressed by following 
expressions : 
x k+1 = F k x k + G k w k 
y k = H k x k + v k 
z k = H k x k 
here , 

x k : a state vector or simply a state, 

w k : a system noise, 

v k : an observation noise, 

y k : an observation signal, 

z k : an output signal, 

F k : dynamics of a system, and 

G k : a drive matrix, 

a maximum energy gain to a filter error from a 
disturbance weighted with the forgetting factor p as an 
evaluation criterion is suppressed to be smaller than a term 
corresponding to a previously given upper limit value y f , and 

the system estimation method, the system estimation 
program for causing a computer to execute respective steps, 
and the computer readable recording medium recording the 
program, includes 
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a step at which a processing section inputs the upper 
limit value y f , the observation signal y k as an input of a 
filter and a value including an observation matrix H k from a 
storage section or an input section, 

a step at which the processing section determines the 
forgetting factor p relevant to the state space model in 
accordance with the upper limit value y f , 

a step at which the processing section reads out an 
initial value or a value including the observation matrix H k at 
a time from the storage section and uses the forgetting factor 
p to execute a hyper Hoo filter expressed by a following 
expression : 

x\ !k = Fk-iX^k-iik-i + K s , k (y k - H k F k _ix" k _i| k -i) 
here , 

x\| k ; an estimated value of a state x k at a time k using 
observation signals y 0 to y k/ 
F k : dynamics of the system, and 
K s , k ; a filter gain, 

a step at which the processing section stores an obtained 
value relating to the hyper Hoo filter into the storage 
section, 

a step at which the processing section calculates an 
existence condition based on the upper limit value yf and the 
forgetting factor p by the observation matrix Hi and the filter 
gain K Sfi , and 

a step at which the processing section sets the upper 
limit value to be small within a range where the existence 
condition is satisfied at each time and stores the value into 
the storage section by decreasing the upper limit value y f and 
repeating the step of executing the hyper H^ filter. 

Besides, according to second solving means of the 
invention, 
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a system estimation device is for making state estimation 
robust and optimizing a forgetting factor p simultaneously in 
an estimation algorithm, in which 

for a state space model expressed by following 
expressions : 
x k+1 = F k x k + G k w k 
y k = H k x k + v k 
z k = H k x k 
here, 

x k : a state vector or simply a state, 

w k : a system noise, 

v k : an observation noise, 

y k : an observation signal, 

z k : an output signal, 

F k : dynamics of a system, and 

G k : a drive matrix, 

a maximum energy gain to a filter error from a 
disturbance weighted with the forgetting factor p as an 
evaluation criterion is suppressed to be smaller than a term 
corresponding to a previously given upper limit value Yf/ 

the system estimation device includes: 

a processing section to execute the estimation algorithm; 

and 

a storage section to which reading and/or writing is 
performed by the processing section and which stores 
respective observed values, set values, and estimated values 
relevant to the state space model, 

the processing section inputs the upper limit value y f , 
the observation signal y k as an input of a filter and a value 
including an observation matrix H k from the storage section or 
an input section, 
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the processing section determines the forgetting factor p 
relevant to the state space model in accordance with the upper 

limit value Yf/ 

the processing section reads out an initial value or a 
value including the observation matrix H k at a time from the 
storage section and uses the forgetting factor p to execute a 
hyper Ho, filter expressed by a following expression: 
x%| k = F k -iX A k -i|)c-i + K s , k (y k - HjcFic-ixVxiic-i) 
here , 

x^ k)k ; an estimated value of a state x k at a time k using 
observation signals y 0 to y k/ 
F k : dynamics of the system, and 
K S/k ; a filter gain, 

the processing section stores an obtained value relating 

to the hyper Hoo filter into the storage section, 

the processing section calculates an existence condition 
based on the upper limit value y f and the forgetting factor p 
by the observation matrix Hi and the filter gain K S(i , and 
the processing part sets the upper limit value to be 
small within a range where the existence condition is 
satisfied at each time and stores the value into the storage 
section by decreasing the upper limit value y f and repeating 
the step of executing the hyper Hoo filter. 

According to the estimation method of the invention, the 
forgetting factor can be optimally determined, and the 
algorithm can stably operate even in the case of single 
precision, and accordingly, high performance can be realized 
at low cost. In general, in a normal civil communication 
equipment, calculation is often performed at single precision 
in view of cost and speed. Thus, as the practical state 
estimation algorithm, the invention would have effects in 
various industrial fields. 
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A more de tailed descr ip tion of the inv en tion is provided 



in the follow ing descript io n and appended claims take n in 

conj u nct ion wi th the accompa ny i ng drawing. 



Brief Dcocr iption of the Drawings brief description of 

THE DRAWINGS 

Fig. 1 is a structural view of hardware of an embodiment. 

Fig. 2 is a flowchart concerning the generation of 
robustness of an Hoo filter and the optimization of a 
forgetting factor p. 

Fig. 3 is a flowchart of an algorithm of the Hoo filter 
(S105) in Fig. 2. 

Fig. 4 is an explanatory view of a square root array 
algorithm of Theorem 2 . 

Fig. 5 is a flowchart of a fast algorithm of Theorem 3, 
which is numerically stable. 

Fig. 6 is a view showing values of an impulse response 

{hih.o 23 . 

Fig. 7 shows an estimation result of the impulse response 
by the fast algorithm of Theorem 3, which is numerically 
stable . 

Fig. 8 is a structural view for system estimation. 
Fig. 9 is an explanatory view of a communication system 
and an echo. 

Fig. 10 is a principle view of an echo canceller. 



DESCRIPTION OF THE PREFERED EMBODIMENTS detailed 

DESCRIPTION OF THE INVENTION 

The following is a detailed description arid explanation 

of the preferred embodiments and best modes for embodying the 
invention along with some proofs thereof. 
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Hereinafter, embodiments of the invention will be described. 

1. Explanation of Symbols 

First, main symbols used in the embodiments of the 
invention and whether they are known or unknown will be 
described. 

x k : State vector or simply a state; unknown and this is an 

object of the estimation. 

x 0 : Initial state; unknown . 

w k : System noise; unknown. 

v k : Observation noise; unknown. 

y k : Observation signal; input of a filter and known. 

z k : Output signal; unknown. 

F k : Dynamics of a system; known. 

G k : Drive matrix; known at the time of execution. 
H k : Observation matrix; known. 

x^ k | k : Estimated value of a state x k at a time k, using 
observation signals y 0 to y k ; given by a filter equation. 
x^ k+ i| k : Estimated value of a state x k+i at a time k+1 using the 
observation symbols y 0 to y k ; given by the filter equation. 
x^ 0 |o: Initial estimated value of a state; originally unknown, 
however, 0 is used for convenience. 

Z^k|k: Corresponds to a covariance matrix of an error of x^ k j k ; 
given by a Riccati equation. 

Z^k+i|k: Corresponds to a covariance matrix of an error of 
x^ k+i | k ; given by the Riccati equation. 

X^ijo: Corresponds to a covariance matrix in an initial state; 
originally unknown, however, e 0 I is used for convenience. 
K s , k : Filter gain; obtained from a matrix X^k|k-i * 
p: Forgetting factor; in the case of Theorems 1 to 3 , when y f 
is determined, it is automatically determined by p = 1 - x(Yf) • 
ef,i: Filter error 
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Re,*: Auxiliary variable 

Incidentally, and M v" placed above the symbol mean 

estimated values. Besides, " - " , "U" and the like are 

symbols added for convenience. Although these symbols are 
placed at the upper right of characters for input convenience, 
as indicated in mathematical expressions, they are the same as 
those placed directly above the characters. Besides, x, w and 
the like are vectors , H, G, K, R, X and the like are matrixes 
and should be expressed by thick letters as indicated in the 
mathematical expressions, however, they are expressed in 
normal letters for input convenience. 

2. Hardware and Program for System Estimation 

The system estimation or the system estimation device and 
system can be provided by a system estimation program for 
causing a computer to execute respective procedures, a 
computer readable recording medium recording the system 
estimation program, a program product including the system 
estimation program and capable of being loaded into an 
internal memory of a computer, a computer, such as a server, 
including the program, and the like. 

Fig. 1 is a structural view of hardware of this 
embodiment . 

This hardware includes a processing section 101 which is 
a central processing unit (CPU), an input section 102, an 
output section 103, a display section 104 and a storage 
section 105. Besides, the processing section 101, the input 
section 102, the output section 103, the display section 104 
and the storage section 105 are connected by suitable 
connection means such as a star or a bus. Known data 
indicated in n l. Explanation of Symbols" and subjected to the 
system estimation are stored in the storage section 105 as the 
need arises. Besides, unknown and known data, calculated data 
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relating to the hyper Hoo filter, and other data are written 
and/or read by the processing section 101 as the need arises. 



3. Hyper Hoo Filter by Which Forgetting Factor Can Be 
Optimally Determined 

(Theorem 1) 

Consideration is given to a state space model as 
indicated by following expressions. 

x k + x = F k x k + G k w k , w k , x k e Tl N (6) 
y k = H k x k + v k , y ki v k en (7) 

z k = H k x kj z k eft, H k G 7l ]xN y * = 0,1,...,L (8) 

An Hoo evaluation criterion as indicated by the following 
expression is proposed for the state space model as described 
above . 

i MVP 

sup £ £ < 7/ 

xo^hw || _ ioi-illV-i + E Kll 2 + E NIVp 

1=0 i=0 ^ 

A state estimated value x^ k)k (or an output estimated value 
z v k(k ) to satisfy this Hoo evaluation criterion is given by a 
hyper Hoo filter of level y f : 
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h\k = H k x k \k 

&k\k = F fc _iXjfc-i|*-i + K Sik {y k - iffcF fc-i&fc— 

= Ek\k-iH' k {HkZ!k\k-\Hl + p)" 1 
E k \ k — £ k \k-i - £ k \ k -iClR~ l k C k X! k \ k -i 



where, 



P 0 



■Re.fc = Rk + C k£ k\k-\C\ , Rk — 

0 < p = l-x(lf) < 1> 7/ > 1 



0 -pyj . 



(10) 

(11) 
(12) 



} (13) 



(14) 
(15) 



Incidentally, expression (11) denotes a filter equation, 
expression (12) denotes a filter gain, and expression (13) 
denotes a Riccati equation. 

Besides, a drive matrix G k is generated as follows. 

GkG l = xh£ Fk £ klkF T (16) 

Besides, in order to improve the tracking capacity of the 
foregoing Hoo filter, the upper limit value y f is set to be as 
small as possible so as to satisfy the following existence 
condition . 

= Et 1 + 1 -^-HjH i >0, i = 0,...,k (17) 

Where, %(yt) is a monotonically damping function of y f/ which 

satisfies = 1 and x(°°) = °- 

The feature of Theorem 1 is that the generation of 
robustness in the state estimation and the optimization of the 
forgetting factor p are simultaneously performed. 

Fig. 2 shows a flowchart concerning the generation of 
robustness of the Hoo filter and the optimization of the 
forgetting factor p. Here, 
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block "EXC > 0" : an existence condition of the Hoc filter, and 

Ay: a positive real number. 

First, the processing section 101 reads out or inputs the 
upper limit value y f from the storage section 105 or the input 
section 102 (S101) . In this example, y f >> 1 is given. The 
processing section 101 determines the forgetting factor p by 
expression (15) (S103) . Thereafter, the processing section 
101 executes the hyper Hoo filter of expression (10) to 
expression (13) based on the forgetting factor p (S105) . The 
processing section 101 calculates expression (17) (or the 
right side (this is made EXC) of after-mentioned expression 
(18) ) (S107) , and when the existence condition is satisfied at 
all times (S109) , y f is decreased by Ay, and the same 
processing is repeated (Sill) . On the other hand, when the 
existence condition is not satisfied at a certain y f (S109) , 
what is obtained by adding Ay to the y f is made the optimal 
value y f op of y f/ and is outputted to the output section 103 
and/or stored into the storage section 105 (S113) . 
Incidentally, in this example, although Ay is added, a 
previously set value other than that may be added. This 
optimization process is called a y-iteration. Incidentally, 
the processing section 101 may store a suitable intermediate 
value and a final value obtained at respective steps, such as 
the Hoo filter calculation step S105 and the existence 
condition calculation step S107, into the storage section 105 
as the need arises, and may read them from the storage section 
105. 

When the hyper Hoo filter satisfies the existence 
condition, the inequality of expression (9) is always 
satisfied. Thus, in the case where the disturbance energy of 
the denominator of expression (9) is limited, the total sum of 
the square estimated error of the numerator of expression (9) 
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becomes bounded, and the estimated error after a certain time 
becomes 0. This means that when y f can be made smaller, the 
estimated value x%j k can quickly follow the change of the state 
x k . 

Here, attention should be given to the fact that the 
algorithm of the hyper H«> filter of Theorem 1 is different 
from that of the normal H*, filter. Besides, when y f -» co, then 
p = 1 and G k = 0, and the algorithm of the H*, filter of Theorem 
1 coincides with the algorithm of the Kalman filter. 

Fig. 3 is a flowchart of the algorithm of the (hyper) Hoo 
filter (S105) in Fig. 2. 

The hyper filtering algorithm can be summarized as 
follows . 

[Step S201] The processing section 101 reads out the initial 
condition of a recursive expression from the storage section 
105 or inputs the initial condition from the input section 
102, and determines it as indicated in the figure. 
Incidentally, L denotes a previously fixed maximum data 
number . 

[Step S203] The processing section 101 compares the time k 
with the maximum data number L. When the time k is larger 
than the maximum data number, the processing section 101 ends 
the processing, and when not larger, advance is made to a next 
step. (If unnecessary, the conditional sentence can be 
removed. Alternatively, restart may be made as the need 
arises . ) 

[Step S205] The processing section 101 calculates a filter 
gain K S/k by using expression (12) . 

[Step S207] The processing section 101 updates the filter 
equation of the hyper Hoo filter of expression (11) . 
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[Step S209] The processing section 101 calculates terms Z"k|k, 
Z^k+ijk corresponding to the covariance matrix of an error by 
using the Riccati equation of expression (13). 

[Step S211] The time k is made to advance (k = k + 1) , return 
is made to step S203, and continuation is made as long as data 
exists . 

Incidentally, the processing section 101 may store a 
suitable intermediate value, a final value, a value of the 
existence condition and the like obtained at the respective 
steps, such as the Hoo filter calculation steps S205 to S209, 
into the storage section 105 as the need arises, or may read 
them from the storage section 105. 

(Scalar Existence Condition) 

The amount of calculation 0(N 2 ) was necessary for the 
judgment of the existence condition of expression (17) . 
However, when the following condition is used, the existence 
of the Hoo filter of Theorem 1, that is, expression (9) can be 
verified by the amount of calculation O(N). 

Corollary 1: Scalar existence condition 

When the following existence condition is used, the 
existence of the hyper Hoo filter can be judged by the amount 
of calculation O(N). 

-gEi + p7/>0, i = 0, (18) 

Here, 

Where, K s ,i denotes the filter gain obtained in expression 
(12) . 

(Proof) 
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Hereinafter, the proof of the system 1 will be described. 
When a characteristic equation 



— H k S k \ k _ x H 1 k A - (— P7/ + H k S k \ k _ l H^) 



\\I-Re*\ = 

= A 2 - (2H k S h]k _ l Hl + p Q )X - pW f + PQHkZ^Hl = 0 

of a 2 x 2 matrix R e ,k is solved, an eigenvalue A, ± of R e ,k is 
obtained as follows. 

. * ± - A P QH k E k \ k .,Hl + 4p2 7 2 
"'Where * = 2^*^-1^ + PQ, Q = 1 " l) 



If 



-4p ff jf *i7 fc | fc _, jf£ + 4/) 2 7 ^ > o 



one of two eigenvalues of the matrix R e ,k becomes positive, the 
other becomes negative, and the matrixes R k and R e , k have the 
same inertia. By this, when 

is used, the existence condition of expression (18) is 
obtained. Here, the amount of calculation of H k K S/k is 0(N). 

4. State Estimation Algorithm Which Is Numerically Stable 

Since the foregoing hyper H*, filter updates Z A k |k-i e R n * n , 
the amount of calculation per unit time step becomes 0(N 2 ), 
that is, an arithmetic operation proportional to N 2 becomes 
necessary. Here, N denotes the dimension of the state vector 
x k . Thus, as the dimension of x k is increased, the calculation 
time required for execution of this filter is rapidly 
increased. Besides, although the error covariance matrix X^kik- 
i must always have the positive definite from its property, 
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there is a case where it has numerically the negative 
definite. Especially, in the case where calculation is made 
at single precision, this tendency becomes remarkable. At 
this time, it is known that the filter becomes unstable. 
Thus, in order to put the algorithm to practical use and to 
reduce the cost, the development of the state estimation 
algorithm which can be operated even at single precision 
(example: 32 bit) is desired. 

Then, next, attention is paid to 



and an H*, filter (square root array algorithm) of Theorem 1, 
which is numerically stabilized, is indicated in Theorem 2. 
Here, although it is assumed that F k = I is established for 
simplification, it can be obtained in the same way also in the 
case of F k * I. Hereinafter, the hyper H*, filter to realize 
the state estimation algorithm which is numerically stable 
will be indicted. 

(Theorem 2) 



Rk-R 1/2 kJ iR^ +/ \> 



R^k— R e,kJlR" 




kjk-l 
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K tJe = K h (:, 1)AM1> 1) . K» = ***** «T Vi«t 



e(*) = 



r i 
Rl 


* I 

C k^kjk-l 


0 


P~* S k\k-l . 





0 




I 



(20) 
(21) 

(22) 



Where, 

i I i 

Rk = -Rj«J i-Rfc . -Rfc = 



r i 



■Re,* = Rk + C r t-^i|t-lCj', C/b = 



0 



1 0 
0 -1 



■27t|t-l - -£*|fc-l-Efc|fc_l 



, R ek = Rl^JiRjj., x 0 |o = Xo 



(23) 



0(k) denotes a J-unitary matrix, that is, satisfies 
0(k)J0H(k) T = J, J = (Ji © I), I denotes a unit matrix, 
K k (:,l) denotes a column vector of a first column of 
the matrix K k . 



Incidentally, in expressions (21) and (22), Ji" 1 and 
Ji can be deleted . 

Fig. 4 is an explanatory view of the square root array 
algorithm of Theorem 2 . This calculation algorithm can be 
used in the calculation (S105) of the Hoo filter in the 
flowchart of Theorem 1 shown in Fig. 2. 

In this estimation algorithm, instead of obtaining £*k|k-i 
by a Riccati type update expression, its factor matrix E^ 1/2 k|k-i 
e R NXN (square root matrix of Z A k|k-i) is obtained by the update 
expression based on the J-unitary transformation. From a 1-1 
block matrix and a 2-1 block matrix generated at this time, 
the filter gain K s , k is obtained as shown in the figure. Thus, 
X^kik-i = I Al/2 k |k-i Z"- i/2 k|k-i > 0 is established, the positive 
definite property of £"k|k-i is ensured, and it can be 
numerically stabilized. Incidentally, a computational 
complexity of the H*, filter of Theorem 2 per unit step remains 
0(N 2 ) . 

Incidentally, in Fig. 4, Ji" 1 can be deleted. 
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First, the processing section 101 reads out terms 
contained in the respective elements of the left-side 
equations of expression (22) from the storage section 105 or 
obtains them from the internal memory or the like, and 
executes the J-unitary transformation (S301) . The processing 
section 101 calculates system gains K k and K g , k from the 
elements of the right-side equations of the obtained 
expression (22) based on expression (21) (S303, S305). The 
processing section 101 calculates the state estimated value 
x^ k | k based on expression (20) (S307) . 

5. Numerically Stable Fast Algorithm for State Estimation 

As described above, a computational complexity of the Hoo 
filter of Theorem 2 per unit step remains 0(N 2 ) . Then, as a 
countermeasure for the complexity, by using that when H k = 
Hk+i^/ H k = [u(k), • • •, u(0), 0, • • •, 0], a covariance matrix 
2 k+ i k | k of one step prediction error of x k = [x T k , 0 T ] T satisfies 



consideration is given to updating L k (that is, L" k ) with a low 
dimension instead of X k+ i| k . Here, when attention is paid to 
R r#k = R 1/2 r , k SR I1/2 r ,k/ next Theorem 3 can be obtained. 

(Theorem 3) 




0 



(24) 
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Xk\k = 2Cjfe-i|fc-i + K a%k (y k - -fir fc £*-ii*-i) (61) 
K,j = K k (: t lyRe^h 1) , K k = pi(K k R;})Rl k (62) 

C k+ iL k R^l 



K k +\ 
0 



e(*) (63) 



here, 0(k) denotes an arbitrary J-unitary matrix, and C k = Ck+i^ is established, where 



1 0 
0 -l 



iie^ = H t + C^^C*', Cfc = J ^* j , K cJk = Rl k J x Rj^ £ 0 (o = *o (23) 



Incidentally, the proof of Theorem 3 will be described 
later . 

The above expression can be arranged with respect to K k 
instead of K" k (= P" 1/2 K k ) . 

Further, when the following J-unitary matrix 

is used, a fast state estimation algorithm of Theorem 4 can be 
obtained. Where, *F denotes a shift matrix. 



(Theorem 4) 
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K Stk = p^K k {: t l)/Re, k {hl) 



K k+ i 




0 


0 







-p *L k R r l k L k C M 



T - T 



0 



ile,*+i = #e,* - Ck+\L k R rk L k C M 

ilr,*+l = -Rr,ik - L k C k ^. x R ek C M L k 



(25) 
(26) 

(27) 

(28) 

(29) 
(30) 



Where, 



ii*+l 



. **+i = («*+i «(* + !- iV)] = K* + l) u*], #i = [ti(l) l 0,... l 0] 



ile.l — -Ri + C'l-S'llO^i > ill = 



1 0 
0 0 
0 1 



10 ~(ry}\ 



, i: 1 | 0 = diag{ / > 2 )/? 3 ,... ) p^}, />=i- x ( 7/ ) 



-1 0 
0 



, = 0, x 0 |o = abo » Kk = p 2 K k (31) 



diagH denotes a diagonal matrix, and R e , k+ i(l,l) denotes a 1-1 
element of the matrix R e# k+i- Besides, the above expression can 
be arranged with respect to K k instead of K~ k . 

In the fast algorithm, since the filter gain K Sfk is 
obtained by the update of lT k e r< n+1 >* 2 i n the following 
factoring 



(32) 



0(N+1) is sufficient for the amount of calculation per unit 
step. Here, attention should be paid to the following 
expression. 







0 


0 







Fig. 5 is one example of a flowchart of a numerically 
stable fast algorithm of Theorem 3-_4. The fast algorithm is 
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incorporated in the calculation step (S105) of the Hoo filter 
of Fig. 2, and is optimized by the y- iteration. Thus, during a 
period in which the existence condition is satisfied, y f is 
gradually decreased, however, at the time point when it comes 
to be unsatisfied, y± is increased as indicated in the figure. 
The Hoo filtering algorithm can be summarized as follows. 

[Step S401] The processing section 101 determines an initial 
condition of the recursive expression as indicated in the 
figure. Incidentally, L denotes a maximum data number. 

[Step S403] The processing section 101 compares the time k 
with the maximum data number L. When the time k is larger 
than the maximum data number, the processing section 101 ends 
the processing, and when not larger, advance is made to a next 
step. (When unnecessary, the conditional sentence can be 
removed. Alternatively, restart is made.) 

[Step S405] The processing section 101 recursively calculates 
a term K k+1 corresponding to a filter gain by using expressions 
(27) and (31) . 

[Step S406] The processing section 101 recursively calculates 
R e ,k+i by using expression (29) . 

[Step S407] The processing section 101 further calculates K s , k 
by using expressions (26) and (31) . 

[Step S409] The processing section 101 judges the existence 
condition EXC > 0 here, and when the existence condition is 
satisfied, advance is made to step S411. 

[Step S413] On the other hand, when the existence condition 
is not satisfied at step S409, the processing section 101 
increases y f , and return is made to step S401. 
[Step S411] The processing section 101 updates the filter 
equation of the filter of expression (25) . 

[Step S415] The processing section 101 recursively calculates 
Rr,k+i by using expression (30) . Besides, the processing 
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section 101 recursively calculates L" k+ i by using expressions 
(28) and (31) . 

[Step S419] The processing section 101 advances the time k (k 
= k+1) , returns to step S403, and continues as long as data 
exists . 

Incidentally, the processing section 101 may store a 
suitable intermediate value and a final value obtained at the 
respective steps, such as the Hoo filter calculation steps S405 
to S415 and the calculation step S409 of the existence 
condition, into the storage section 105 as the need arises, 
and may read them from the storage section 105. 

6. Echo Canceller 

Next, a mathematical model of an echo canceling problem 
is generated. 

First, when consideration is given to the fact that a 
received signal {uk} becomes an input signal to an echo path, 
by a (time-varying) impulse response {hi[k]} of the echo path, 
an observed value {y k } of an echo {d k } is expressed by the 
following expression . 

Vk = d k + v k = *52hi[k]u k -i+v k , fc = 0,l,2, • • (33) 

Here, u k and y k respectively denote the received signal and the 
echo at a time t k (= kT; T is a sampling period) , v k denotes 
circuit noise having a mean value of 0 at the time tk, h± [k] , i 
= 0, * • •, N-l denotes a time-varying impulse response, and the 
tap number N thereof is known. At this time, when estimated 
values {h\[k]} of the impulse response are obtained every 
moment, a pseudo-echo as indicated below can be generated by 
using that. 
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/v-i 

dk = 52hi[k]uk-i, fc = 0, 1, 2, • • (34) 

When this is subtracted from the echo (y k - d^ k » 0), the echo 
can be cancelled. Where, it is assumed that if k-i < 0, then 
u k -i = 0 . 

From the above, the problem can be reduced to the problem 
of successively estimating the impulse response {h±[k]} of the 
echo path from the received signal {u k } and the echo {y k } which 
can be directly observed. 

In general, in order to apply the Hoo filter to the echo 
canceller, first, expression (32) must be expressed by a state 
space model including a state equation and an observation 
equation. Then, since the problem is to estimate the impulse 
response {hi[k]}, when {hjjkj} is made a state variable x k , and 
a variance of about w k is allowed, the following state space 
model can be established for the echo path. 

x fc+ i = x k + GkW kl x k ,w k en N (35) 

Vk = H k x k + v k , yk,v k en (36) 

z k = H k x k , z k en, H k eK lx » (37) 

Where, 

x k = [M*]i---. fc ^i[*]] r , w k = [w k (l)r-<,Wk(N)} T 
H k = [tijb, • • • ,ujb-iv+i] 

The hyper and fast filtering algorithms to the state 
space model as stated above is as described before. Besides, 
at the estimation of the impulse response, when the generation 
of a transmission signal is detected, the estimation is 
generally suspended during that. 

7. Evaluation to Impulse Response 

(Confirmation of Operation) 
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With respect to the case where the impulse response of 
the echo pulse is temporally invariable (hi [k] = hi), and the 
tap number N is 48, the operation of the fast algorithm is 
confirmed by using a simulation. 

* = + (38) 

Incidentally, Fig. 6 is a view showing values of the impulse 
response {hi} here. 

Here, the value shown in the figure are used for the 
impulse response {hi} i=0 23 , and the other {hi} i=2 4 47 is made 0. 
Besides, it is assumed that u k is stationary Gaussian white 
noise having a mean value of 0 and variance a v 2 = 1.0 x 10" 6 , 
and the sampling period T is made 1.0 for convenience. 

Besides, the received signal {u k } is approximated by a 
secondary AR model as follows. 
u k = aiu k _i + a 2 u k _ 2 + w k * (39) 

Where, oti = 0.7 and a 2 = 0.1 are assumed, and w k ' denotes 
stationary Gaussian white noise having a means value of 0 and 
variance a w . 2 = 0.04. 

(Estimation Result of Impulse Response) 

Fig. 7 shows an estimation result of the impulse response 
by the numerically stable fast algorithm of Theorem 43-. Here, 
the vertical axis of Fig. 7(b) indicates 
^{Zi^Chi-x^Ci + l)) 2 }. 

By this, it is understood that the estimation can be 
excellently performed by the fast algorithm. Where, 
P = 1 - X(Yf)/ X(Yf) = Yf" 2 ' x ^o|o = 0 and Z'ojo = 201 were assumed, 
and the calculation was performed at double precision. 
Besides, while the existence condition is confirmed, y f = 5.5 
was set . 
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8 . Proof of Theorem 



8-1. Proof of Theorem 2 

When the following expression: 



0 p~*El\k-i j 



T 

JtJ 0 



Rl 



0 



4 — l "5 



j4 p-^^xt 

0 



(40) 



is established, following expressions are obtained by 
comparing the respective terms of 2 x 2 block matrixes of both 
sides . 

R ek = Rk + Ck£k\k-\C£ (41) 
K k = iJ klk - X Cl (42) 
+ p~ l K f^R'^K^ = p" 1 !?^^ (43) 

This is coincident with the Riccati equation of 
expression (13) at F k = I of Theorem 1. Where, 



J= (Ji0/), J, 



-[:-,]■ 



C k = 



(44) 



On the other hand, when AJA T = BJB T is established, B can 
be expressed as B = A0(k) by using the J-unitary matrix 0(K) . 
Thus, from expression (4 0) , the Riccati equation of Theorem 1 
is equivalent to the following expression. 



o 



2 



R k C k El\ k _ x 



e(*) 



(45) 



Incidentally, in expressions (40) and (45) , Ji" 1 can be 
deleted. 



8-2. Proof of Theorem 3 
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It is assumed that there is a J-unitary matrix 0(k) which 
performs block triangulation as follows. 



At this time, when both sides J = (Ji©-S) - norm of the above 
expression are compared, X, Y and Z of the left side can be 
determined as follows. Where, S denote a diagonal matrix in 
which diagonal elements take 1 or -1. 
(1,1) -Block matrix 



X 0 
Y Z 




XJ X X T 



= R** + c k+l (i? t+1 , fc - *i7 t | t _ 1 * r ) 



• * = is obtained from = R^J^^, = 



Here, attention should be paid to the fact that 



j r i » j^ji = d, s- 1 = s, - = Rrt, = r;*sr;}, 



is established. 



(2,1) -Block matrix 
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YJ,X T 



k+l 



By this, Y _ ^ -Kfc+i | jRj* +1 J, is obtained. 
Where, cl = (C k+l #) T 



(2, 2) -Block Matrix 
-ZSZ T + YJ 1 Y T 

= ([ ] R$ [%']'- i^) ** + p->£ k+l[k 

By this, 

and Z = i Jfe+1 H r r^ 1 is obtained. 



8-3. Proof of Theorem 4 

When an observation matrix H k has a shift characteristic 

and 

J = ( Ji © - S) , 

the following relational expression is obtained by a similar 
method to Theorem 2 . 
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Where, 



K k +\ 

0 



0 



0 



E<*> 



e(*)= (j.rI e-/£)£ MiKlS? ®-<* 2 *.) 



J7(*) = 



Where, L(fc) r (R Ciit © K*) = ® -«r.*+i) 



(46) 



(47) 



and R r ,k+i is determined so that Z(k) T (R e ,k © - Rr,k)Z(k) = 
(R e ,k+i ® ~ Rr,k+i) is established. Next, when an update 
expression of R r ,k+i is newly added to the third line of 
expression (46), the following expression is finally obtained. 



K k +\ 
0 
0 



Where, 



0 

Rr,k+\ 



Ck+\L k 

Rrjt 



Re,k 

0 

- T « T 

L k Cfc+i 

Re,k ~ C k +iL k R~lL k C 
0 

Kk 



-R-elCk+lLk 



fc+1 
T - T 



-p ^L k R rJc L k C k+l p *L k - 



K k 



R ei lCk+iL k 



Rr,k - L k C k + x R ek C k +\L 



(48) 



From the correspondence of the respective terms of 3 x 2 
block matrixes of both sides, the following update expression 
of a gain matrix K" k is obtained. 
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0 


0 







p -tL k R;lL' k C' k+l 



0 



Rr^+i = Rr,h ~ L k C k+1 R~ l k C k +iL k 



(49) 

(50) 

(51) 
(52) 



Industrial Applicability 



In general, in a normal civil communication equipment or 
the like, calculation is often performed at single precision 
in view of the cost and speed. Thus, as the practical state 
estimation algorithm, the present invention would have effects 
in various industrial fields. Besides, the invention can be 
applied to an echo canceller in a communication system or a 
sound system, sound field reproduction, noise control and the 
like. 

Although embodiments of the invention have been shown and 
described, it is to be understood that various modifications 
and substitutions, as well as rearrangements of method steps 
and equipment, parts and components, can be made by those 
skilled in the art without departing from the novel spirit and 
scope of the invention 
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